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1. Abstract 



A novel computational treatment of dense, stiff, coupled reaction rate equations is introduced 
to study the nucleation, growth, and possible coalescence of cavities during neutron irradiation 
of metals. Radiation damage is modeled by the creation of Frenkel pair defects and helium im- 
purity atoms. A multi-dimensional cluster size distribution function allows independent evolu- 
tion of the vacancy and helium content of cavities, distinguishing voids and bubbles. A model 
with sessile cavities and no cluster-cluster coalescence can result in a bimodal final cavity size 
distribution with coexistence of small, high-pressure bubbles and large, low-pressure voids. A 
model that includes unhindered cavity diffusion and coalescence ultimately removes the small 
helium bubbles from the system, leaving only large voids. The terminal void density is also re- 
duced and the incubation period and terminal swelling rate can be greatly altered by cavity co- 
alescence. Temperature-dependent trapping of voids/bubbles by precipitates and alterations 
in void surface diffusion from adsorbed impurities and internal gas pressure may give rise to 
intermediate swelling behavior through their effects on cavity mobility and coalescence. 



2. Introduction 



Irradiation of metals has long been known to culminate in macroscopic property changes in- 
cluding void swelling pQ. Characteristic stable voids and steady volumetric swelling develop 
for a range of temperatures and fluxes, independent of whether radiation bombardment dam- 
age occurs as disseminated Frenkel pairs or as small defect clusters. This can occur whether or 
not helium is generated along with atomic displacements. In either case, small, unstable voids, 
loops, and other defect clusters will develop almost immediately within the irradiated material. 
Their subsequent evolution determines the fluence required to create stable voids and achieve 
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steady swelling; this so-called incubation dose includes most of the dependence on radiation 
environment [2f3"P] . The processes that govern microstructure evolution include thermally- 
activated motion of small defect clusters, mutual impingement, and annihilation or coalescence 
reactions along with micro-chemical changes from nuclear transmutation and displacements 
or diffusion of pre-existing impurities. Radiation simulations should ideally encompass all of 
these processes. Typically, existing models have included only particular types of defects and 
reactions or have made other numerical approximations in order to obtain a solution. 

At the least, simulations of early irradiation must account for void nucleation and growth pro- 
cesses, since annihilation, aggregation, and cluster ripening take place concurrently. Transient 
and steady-state swelling behavior due to these processes have been studied recently [5116117118] . 
However, only void reactions with vacancy or interstitial monomers are included in these stud- 
ies. This minimal model of void nucleation gives reasonable swelling behavior as a function of 
temperature and flux [7f5] . viz. an observed steady swelling rate around 1%/dpa in austenitic 
stainless steels and an important flux-effect on the measured incubation times pfTU] . While the 
results are encouraging, these calculations neglect many of the processes believed to shape the 
microstructure. For example, the generation and aggregation of helium impurities is not ex- 
plicitly modeled. Size-dependent void diffusion [111112] is neglected, and thus direct void-void 
coalescence is not included. Dislocation loop formation, migration, and coalescence is not ex- 
plicitly simulated, either. The model can be considered to combine the production and biased 
diffusion of small vacancy and interstitial clusters into effective generation and reaction rates 
for monomer species alone, but it is unclear a priori how a coarse-grained treatment of these 
processes affects microstructure evolution. 

It is now clear that the model must presuppose a ready supply of gas impurity atoms (e.g., 
oxygen and helium [13J) to promote the formation of small voids from the radiation-induced, 
supersaturated vacancy population. In practice, reasonable corrections to void energies may 
reproduce the approximate void number density observed in irradiated steel [8]. Ultimately, 
however, crude models for the apportionment of impurities among the defect clusters should be 
supplanted by a detailed accounting of multicomponent aggregation and coalescence reactions 
and their influence on the non-equilibrium cluster size distribution. Such problems are widely 
addressed in the literature, including gelation, polymerization, and formation of aerosols and 

numerical methods developed for such problems may also be fruitfully applied to radiation 
swelling. Here, a hybrid numerical approach that can encompass the full range of possible clus- 
ter compositions and cluster reactions in mean field is introduced, a Livermore Microstruc- 
ture Evolution program, LiME. As a first application, the method is applied to the nucleation 
and growth of voids with a two-component distribution of cluster compositions, examining the 
evolution of helium- vacancy clusters |13j . while continuing to treat oxygen adsorption by sim- 
ply reducing the cavity surface energy by a constant (temperature-independent) factor. The 
method predicts realistic swelling behavior for ferritic steel in reactor environments. 

As before, the void distribution function is partitioned into overlapping regions [5], treating 
small clusters with the Master Equation (ME domain) and large ones with Monte Carlo meth- 
ods (MC domain). This allows self-consistent evolution of the full void population with no 
truncation of the size domain, no assumptions as to the critical void size or the nature of the 
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nucleation process, and no approximations for the overall nucleation rate or duration of the 
nucleation stage. Monomer concentrations are included in the ME region, where they may ei- 
ther be treated separately by a quasi-stationary approximation or evolved along with the small 
clusters through coupled nonlinear reaction rate equations. The formation and evolution of dis- 
location loops is not explicitly modeled; network dislocations and loops are already described 
by a single, time-dependent density parameter rather than a detailed size distribution function 
[3T] . However, the methods used for void evolution would be easily generalizable to other defect 
species and reactions, provided that suitable mean field rate coefficients are specified for their 
reaction rate equations. In particular, future calculations will consider the formation, unfault- 
ing, and migration of dislocation loops; loop coalescence and annihilation; and incorporation 
of loops in the dislocation network. 

The remainder of this paper first describes the coupled, stiff, non-linear evolution equations for 
void nucleation, growth, and coalescence. It presents the microscopic rate theory model, gives 
an overview of the computational scheme, details the various numerical methods employed in 
the calculations, and makes a preliminary application to void nucleation in irradiated stainless 
steel. The simulations include vacancy, interstitial, and helium generation, aggregation and, 
annihilation, with or without cluster coalescence. The results are sensitive to the effects of 
absorbed impurity atoms on cavity surface energy. They also expose a substantial influence of 
small, unstable, mobile clusters on the formation of critical-sized voids via direct cluster-cluster 
coalescence. Realistic incubation and swelling behavior cannot be obtained over wide ranges 
of temperature and flux without including cluster mobility and coalescence. 



3. Rate Theory Model 

Allowable microstructure reactions (either aggregation or annihilation) are assumed to occur 
whenever two defects, m and n, come into contact. Within the mean field continuum approxi- 
mation, the collision rate is proportional to their relative diffusivity, D m n , and effective colli- 
sion cross-section, A m n . As before [3]f5] . a bias factor Z m n includes the effect of long-range in- 
teractions [321(33] on the binary reaction rates, K(m, n)p m p n , where p are densities of reactant 
species m 7^ n and the rate coefficients are: 



Note that an additional factor of 1/2 may be required when m = n, to prevent double-counting 
of unique pairs of identical reactant particles. This factor is not explicitly shown in the definition 
ofK. 

Microstructure defect species are limited here to self-interstitials and -vacancies, substitu- 
tional and interstitial helium, voids/bubbles, and network dislocations. Vacancy and helium 
monomers as well as clusters are characterized by their composition, n = (n vac , Uhei), in a two 
dimensional space. Self-vacancies and interstitials are also specially identified by (1, 0) = v and 
(—1,0) = i, respectively; substitutional and interstitial helium by (1, 1) = vh and (0, 1) = h; 
and network dislocations by d. Monomer densities evolve according to: 
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The vacancy- vacancy aggregation term, (Z VjV A V)V D VjV )p v p v , within the first summation for 
dp v / dt in Eq. [2]includes that two vacancies are consumed by the reaction, that there is a factor 
of 1/2 to prevent double-counting of unique pairs of vacancies from the population p v , and that 
the relative diffusivity is twice D v . The net rate is identical to that used in a previous study [5] . 
Similar considerations also apply to pairs of substitutional helium and to thermal dissociation 
of vacancy dimers. 

Cluster (n ^ {v, vh, h, i}) densities evolve as: 
dp n 



dt 
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(3) 



in terms of any direct generation of clusters in the radiation damage cascade, g n ; reactions of 
existing clusters with monomers (in brackets) that consume or create n-mers including thermal 
emission of vacancies, and cluster-cluster reactions (in the second set of brackets) that consume 
or create n-mers. Factors of 1/2 in the first and last summations prevent double counting of 
indistinguishable reactant pairs, and <5 mjI1 = S mvt n v 5 mhtnh where the right hand side consists 

!1 i — j 
. The primed summation is restricted to all 
im- 
pairs of reactants with m, n — m ^ {v, vh, h, i}. Defects n — m (n — v , etc.) are restricted to 
the domain of allowed compositions by a step function: U(n) = U{n v )U{n h ), where U{n) = 

1 n > 

~~ . Finally, clusters never undergo fission in this model, only the thermal emission of 
n < 
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single vacancies. 



Radiation damage deposition is approximated by the creation of disseminated monomers, so 
g n = for n ^ {v, vh, h, i}. In this case, g^ = 0£, in terms of the atomic displacement rate, 0, 
and the damage production efficiency, £. The total helium production is g^ + g v h, with the ratio 
of interstitial to substitutional depending on the model. (Here, it is assumed that the helium is 
all deposited as substitutional defects.) Conservation of host atoms (including transmutation 
products) requires g v +9vh = 9%- Helium impurities are added with a temperature-independent, 
gradual activation of a-emitters. This is modeled for a Fe-Ni-Cr steel undergoing neutron bom- 
bardment according a two-step activation process, in analogy to the 58 Ni(n,7) 59 Ni(n,a) reac- 
tion. Model transmutation rates are treated as free parameters and are fit to the experimental 
helium content in HFIR-irradiated nickel p4T3"5] . The parameters are 7, a, and 5 for the rates 
of (respectively) 58 Ni(n,7), 59 Ni(n,a), and the sum of all transmutations that consume 59 Ni. In 
terms of the cumulative radiation dose in dpa, x — J <p(f)dt (for radiation flux, <p): 



d 

dx 



' ' P58 
U 9 . 




^58 
U 9 , 



(4) 



The 59 Ni content, p^g, is obtained from Eq.[4]by transforming to the eigenbasis, where Pa{%) = 
P5s{x) and Pb(x) = p^(x) + ^psPssix) are solved, and then transforming back. The helium 
generation rate is given by: 

^ = ap 59 (x) = a^-p 58 (0) \e- Sx - e~H (5) 
dx 7 — 1 J 

assuming that only 59 Ni(n,a) produces a-particles. The fit parameters are 7 = 0.0255, a = 
0.0711, and 6 = 0.297 dpa -1 . Pristine type 316 stainless steel is approximately 14% nickel, with 
68.08% of that 58 Ni and with no naturally-occurring 59 Ni. Other relevant materials parameters 
for type-316 stainless steel are listed in Table [Tj 

Non-interacting diffusion (independent random walks) implies D m n = D m + D n . Defect col- 
lision cross-sections are simply given by 

A m>n = 47r(r m + r n ) for m ^ {v, vh, h, i} and n ^ {v, vh, h, i} 

A m ,n = r m + b for n G {v, vh, h, i} (6) 



in terms of radii for (spherical) defects, r n = y ( n ^ ' (except for interstitial monomers, where 
Ti = Th = r v ). For consistency with earlier work, cross-sections involving monomers are defined 
using the Burgers vector magnitude in place of a monomer radius. 

Bias factors between voids and the four defect monomers are calculated from a mean field so- 
lution of the diffusion including stress-mediated interactions (33] . The infinite series describing 
the image interaction [36] is fit by a simple analytic form, while the modulus interaction [37] 
is treated analytically. The numerical results are tabulated for small voids and computed as 
needed for larger ones. Long range void-void interactions are presently neglected, so Z m n = 1 
for m, n ^ {v, vh, h, i}. In principle, the effect of any long-range interactions or net drift ve- 
locities (e.g., from external stress or temperature gradients [38J) can be incorporated in the 
void-void reaction rates, so the mean field reaction kernel, K, has general applicability. 
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Thermal emission from vacancy clusters is evaluated by a detailed balance argument. Equating 
vacancy emission and absorption for the n-mer identifies the chemical potential, /j,]^ = — 
p[n-v]^ j n t erms of the n-mer and (n — t>)-mer (i.e., void minus one vacancy) free energies. 
Rewriting in terms of void internal energies, E, and the inert gas pressure, P: 

C>] = c [eq] e (E^-E^-Pn)/kT ^ 

Gas pressure is described with a non-ideal equation of state for helium versus density and 
temperature [39] . No volume relaxation is included (i.e. , the void volume is n v Q) . In the absence 
of surface-adsorbed impurity atoms, E^ is parametrized in terms of an effective surface energy, 
7'"', and the surface area of a spherical cavity of volume n v VL 

E [n] = 7 H 47rr 2 = A 7 0( T ) U _ 4yrr 2 (g) 

In the continuum limit, 7^ approaches that of a flat, clean surface, 70(7"), while it approaches 
the results of atomic calculations in the limit of small voids [30] . This surface energy is then 
further reduced by an constant scale factor, A, to reflect the presence of adsorbed oxygen 
impurities jH] (see Table [l|. Finally, the emission rate is obtained from c[ n ' and the vacancy- 
cluster reaction parameters for the (n— v)-mer. For straight, jogged dislocation segments, cjfl = 
cjf^, the thermal equilibrium concentration. Emission rate coefficients in Eq.|2]are represented 
as unary reactions, by including the defect-dependent c{p within the rate coefficient. 

At some maximum density, an over-pressurized bubble would begin to emit self-interstitials via 
loop punching [39] • Such a possibility is not considered here; instead, an artificial constraint is 
imposed on the helium density in a reactant cluster, < 2n v . Any reactions that would yield a 
higher density are disallowed. Thermal dissociation of substitutional helium into a self- vacancy 
plus interstitial helium is similarly assumed to be energetically impossible at temperatures of 
interest. Note that self-interstitial and interstitial helium aggregation is excluded since inter- 
stitial loops are effectively part of the dislocation density model. Mixed interstitial clusters can 
develop in principle [32] ■ 

Void diffusivity D n = D v /n v i ^ 3 for n = {n v ,rih). This gives both the correct monomer value 
and size-dependence for large cluster diffusion [TT1IT2"] . although the activation energy for void 
migration should more properly be that for surface diffusion. This diffusivity takes no account of 
the effect of reversible pinning [13], or internal gas pressure on the migration [33], or radiation- 
enhanced diffusion [35], or, e.g., that vacancy dimer diffusion may be D V2 ~ D v . Trapping at 
dislocations and grain boundaries are not considered. Such features would be straightforward 
to incorporate in the future. 

The dislocation model reproduces measured dislocation densities versus dose and temperature 
|31j . It includes separate source and annealing terms in terms of the biased flow of radiation- 
induced vacancies and interstitials. There is one adjustable parameter, /, representing a char- 
acteristic dislocation pinning length [31]. This is taken to be independent of the density of 
voids/bubbles in the matrix, because the pinning length in stainless steels is more determined 
by carbide nano-precipitates than by the density of voids/bubbles. 
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4. Numerical Method 



4.1. Overview 



Once the temperature- and radiation-environment are specified and initial conditions for the 
microstructure are fixed, the Master Equations[2]and[3]completely determine the void/bubble 
size distribution function P(t) = {p n (t)}. Such stiff, non-linear coupled rate equations can be 
integrated numerically [TT] , although this becomes intractable for a large domain of cluster 
sizes. The number of distinct species may be reduced by grouping similar clusters together [46J , 
but the direct approach still becomes intractable for multi-dimensional distributions. Monte 
Carlo schemes for discrete coalescence events [T5|20] can naturally encompass large voids of 
arbitrary composition; however, they are inefficient for simulating nucleation from sub-critical 
clusters. Here, the advantages of both methods are combined by partitioning the cluster compo- 
sition domain into two overlapping regions. Separate sub-distributions are defined for each, la- 
beled ME and MC for treatment by Master Equation and Monte Carlo, with P = P ME + P MC . 
Each sub-distribution is composed of discrete ensembles of identical clusters, represented by 
(n, p) for the paired multi-dimensional cluster composition, n, and the ensemble density, p. 
The distribution p ME = {(n, p)} ME includes interstitials, i, and vacancy-helium clusters, n, 
with < n v < N™ E and < n h < Nff E . There is exactly one element for each ME species, 
for a total of N ME . Only the densities of the ME elements evolve over time. A sparse, random 
set, {(n, p)} , approximates P MC for all < n v ,nh < oo. The total number of elements, 
N MC , is variable, and there may be none, one, or many MC elements for a given n (each with 
potentially different values of p). Both the densities and the compositions of the MC elements 
evolve with time. Such split distribution functions have been used before in a Fokker-Planck 
treatment of void growth [5], and in non-equilibrium chemistry [23.29], and plasma physics 
applications [47] . In essence, the elements of P MC also constitute so-called "macroparticles", 
already in wide use for non-equilibrium plasma physics problems [48J . 



ME-ME reactions (those processes with reactants and product among the elements of P ME ^ are 
evaluated in a continuum approximation, using the Master Equation [TTJ. Discrete MC-MC re- 
actions are performed stochastically using a Markov Monte Carlo procedure [TS] . ME-MC cross- 
reactions are included using either the Markov Monte Carlo method or Poisson-distributed 
random walks [2T)]|2~1~] for p MC ; and using average sink or source terms in the rate equations 
for p ME . There are also procedures to transfer clusters between the two sub-distributions and 
to regulate the number of elements and their ensemble densities in P MC ; in order to control 
statistical errors and computational cost. This mixed algorithm is elaborate, so the different 
approaches for each of the various components are described in detail in the following sections. 

The material microstructure is evolved over time-step, r, by operator splitting into five stages. 
First, ME-MC reactions for rapidly evolving MC clusters (i.e., those with small n v ) are included 



(Sec. 4.5) with a Markov chain method. Second, the ME-MC reactions for the large, slowly- 



evolving clusters are evaluated by Poisson-distributed random walks in composition space for 



each possible reaction with ME species (Sec. 4.5). Third, all MC-MC reactions are evaluated 
with the Markov Monte Carlo method. (Sec. 4.4) This completes the evolution of P MC over 
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r. The fourth stage integrates the ME including the average source and sink terms from MC 



defects and dislocations (Sec 4.2 ). This completes the evolution for the void/bubble P. At this 



point, clusters may be exchanged between p ME and P MC \ without affecting the instantaneous 



total P in any way (Sec. 4.3 ) . This procedure may create new MC elements or eliminate existing 
ones, in order to control the growth of N MC versus time. Fifth and finally, dislocation evolution 
is performed using a previously-described model [31] . 



Overall numerical accuracy is monitored through the conservation of host and helium atoms. 
Operator splitting of the evolution equations causes first-order time-step errors. However, con- 
servation errors are dominated by differences between the ME and MC treatment of reac- 



tions in Sec. 4.5 (i.e., continuous reactions versus discrete, stochastic events). These artifact 
statistical fluctuations are most important at low temperatures and especially during incuba- 
tion, when N MC is smaller, defect annihilation dominates, and little net swelling occurs. They 
must be carefully controlled, since the transient period represents a sort of barrier-crossing 
problem, with nucleation of stable voids and concomitant, self-consistent changes in the va- 
cancy/interstitial populations as the barrier. Any artificial Monte Carlo noise must not spuri- 
ously affect the crossing into the steady-state. In other words, -/V MC, -dependent fluctuations in 
the net vacancy content must not significantly promote or inhibit void nucleation. In practice, 
stable cavities form naturally under the vacancy supersaturation and essentially irreversible 
aggregation of helium, and the volumetric swelling behavior is not unduly sensitive to N MC 
for the situations considered here. 



4.2. ME-ME reactions 



Small defect clusters develop at high concentrations under irradiation, and so they dominate 
the system of reactions. However, they quickly reach a quasi-stationary distribution wherein 
further reactions cause little change in their densities; i.e., the majority of their reactions sub- 
sequently cancel one another. It is much more efficient to treat the net reaction rates in a con- 
tinuum approximation rather than to explicitly account for individual reactions. The ordinary 
differential equation solver, VODE, provides an optimized treatment of stiff, nonlinear reaction 
equations [T7] , given f n =^f (Eqs. 2 and 3) and the Jacobian, J nm = J^ 2 - for all species. The 
computational cost increases rapidly with the number of coupled equations, hence the cluster 
domain is limited to < n vac < N™J and < n hel < N^f. Typically, N**J? = 10-100 and 
Nfai = 2-10. Some terms are excluded from the Master Equation so that all reaction products 
remain within this finite domain. Clusters with < n vac < iV^f /2 and < n^i < Nj^f /2 
may undergo any mutual reactions, but no other ME clusters may undergo any reactions. These 
latter clusters are frozen in size, so their density only increases as reaction products accumu- 



late. Frozen clusters eventually transfer to the MC distribution as described in Section 4.3 
after which they will undergo reactions normally. 

With reaction constraints and separate ME and MC distributions, the vacancy Eq.[2]becomes: 
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dfv 
dt 



g v (t) + K(v 2 ,i)p V2 {t)pi{t) 

+ ]T -K(n,v) Pa (t)p v (t)+K(n,0)p n (t) 



U(-N ME 



n 



S fast + S slo W{to) pv{t) + gfast + gslon,^ 



(9) 



restricting the sums over n G ME to reactive defects. Eq. [TJalso parametrizes unary vacancy 
emission reactions as the n-null reaction, K(n,0) = Z n _ v v A n _ VjV D n _ VjV c)* l ~ v \ S includes the 
external source and sink terms for reactive elements of P ME ; it accounts for ME reactions 
with defects in p MC and with dislocations. Vacancy absorption at MC defects and dislocations 
is parametrized by S v , and vacancy emission by So. The vacancy sinks and sources include 
separate terms that either evolve slowly or rapidly with time. The coefficients are obtained in 



Sec. 4.5 The rest of Eq. Stakes similar form, with sinks Si, S v h, or Sh- (Only vacancies can be 



thermally emitted from defect clusters, so So is the only source term.) 

Operator splitting over the time-step, r, is such that external source and sink terms S are held 
constant as P ME evolves. S is divided into terms that evolve slowly or rapidly with time. The 
bar indicates an average of the sink strength over the time-step, from t to to + T i useful for 
rapidly evolving MC clusters, while slowly-evolving dislocations and large MC voids are simply 



evaluated at the beginning to (see also Sec. 4.5 for further details). 
The constrained coalescence Eq. [3]becomes: 
dp n { ^ 5 r 



dt 



9n(t) + < Y K(m,n - m)/o m (t)p n _ m (t)(l 
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(10) 
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for clusters m, n G ME, and n {v,vh, h,i}. The primed summation excludes n — m G 
{v, vh, h, i}, since the monomer reactions are evaluated separately. S includes any reactions of 
the n-mer with the MC clusters and with dislocations. There are no reactions that consume 
frozen clusters, so their concentration increases with time. 



A subset of the disallowed reactions would produce clusters that still lie within the ME domain. 
These have been excluded, for simplicity and to better resemble an earlier scheme for monomer 
aggregation [5j. Specifically, a homogeneous boundary condition is imposed on the Fokker- 
Planck equation in Ref. [5], at n = N^f. Clusters that grow to the boundary are removed 
from the Master Equation treatment and accumulated separately, during which time they are 
prevented from changing size. This is equivalent to keeping those iV^f-sized clusters within 
P ME but disabling all of their reactions. Frozen clusters are then intermittently transferred to 
P MC> , where they are no longer constrained [SJ. 

Ideally, the ME domain will encompass all non-zero generation terms, g n , and include as many 
sub-critical or transient defect cluster species as possible. A relatively small domain of N^f E ~ 
60, Nff E ~ 4 is chosen here, reflecting the computational demands that coalescence imposes. 
Similarly to |5J, the solution is recorded at exponentially-increasing intervals. This time-step 
is irrelevant to the ME evolution itself, which advances by adaptive sub-steps. However, r 
controls errors from operator splitting of the evolution equations, and it governs the creation 
of MC elements, as described below. 

Because the sink/source terms, S, are evaluated by a discrete MC method, they introduce a fic- 
titious noise to the continuum rate equations. This partly manifests as step-function disconti- 
nuities in the sink strength over successive time-steps, which in turn causes transient relaxation 
in the concentrations of the ME species. The numerical solution tries to accurately follow the 
transients, potentially making the fully coupled, non-linear evolution inefficient, when large 
time-steps are otherwise possible. Rather than faithfully simulating these spurious transients 
at late times, it may be preferable to solve the monomer concentrations (Eqs. |2j [9j etc.) in the 
quasi-stationary approximation after any real transient behavior (due to the abrupt onset of 



irradiation or other changes in environmental parameters) has concluded. Eq. 10 for dimers 
and larger clusters may then be solved while holding the monomer concentrations fixed over 
the time-step. In practice, after a brief transient, the results are comparable to those obtained 
from the full, coupled, non-linear ME solution. 



4.3. Transfer between ME and MC domains 



A majority of the ME elements in a small multi-dimensional domain will lie near its boundary, 
and so the majority of the ME cluster species will be artificially frozen. The constraints on the 
defect clusters are only lifted after they are transferred to P MC . There are three desiderata to 
this transfer process. Foremost, it must minimize any systematic, constraint-induced errors, 
therefore the density of frozen clusters must be small compared to the rest of P. Secondly, the 
MC computational cost must be controlled, therefore N MC must be kept small. Rather than 
increasing N at every opportunity, frozen clusters at n e ME are allowed to accumulate 
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until exceeding a spawning threshold density, p^ IE > p sp , as in [5] . At the end of that time-step, 
a portion of the accumulated density is removed from p ME and transferred to a new element 
of P MC , incrementing N MC . 

(n, p n ) ME -> (n, p n - 6p) ME + (n; 6p) MC (11) 

with the ME and MC compositions coinciding. If the accumulated p^ IE > p sp after each time- 
step, then the accumulating clusters are effectively never constrained. Finally, it is imperative 
to minimize any iV A/c -dependent Monte Carlo statistical error. Individual MC elements with 
the largest p will contribute the most to this error. Therefore, if p^ E 3> p S p at the end of a 
time-step, then AN > 1 new MC elements are created, as: 

(n, p n ) ME - (n, Pn - 5p) ME + AN x (n; ^) MC (12) 

Equivalently, MC elements with large p may be split into identical macroparticles with smaller 
densities. The chosen values for r, p sp , and the functional dependence of AiV on 5p control the 
N -related statistical error and computational cost for a simulation. Typically, log 2 (AiV) = 
Int(log 30 ((5p/p S p)). 

For example, the distribution in Fig. [TJshows the production of many MC macroparticles con- 
taining 2-4 helium atoms; these react and form a plume that extends to n v ~ 100. The ME 
domain used in this example also includes frozen cluster species with 5-9 helium; these species 
have not yet reached the threshold density. They eventually spawn MC elements, but at a much 
slower rate than for the near-critical sizes of 2-4 helium. Even at this early time, the total den- 
sity of constrained ME clusters is small compared to the MC population so constraint errors 
are minimized. 



Since p sp cannot be made arbitrarily small in practice, it is useful to add a second transfer 
mechanism. When a pre-existing MC element at n falls inside the frozen ME domain, the 
change: 

(n, p n ) ME + (n, p n ) MC -> (n, p n - Sp) ME + (n, p n + Sp) MC (13) 

leaves the total distribution unchanged. N MC remains constant, so the calculation remains 
tractable. In practice, the maximum amount 5p < p n is transferred until the receiving MC 
element reaches a cutoff density, p^ IC + 5p < p max (where typically, p max ~ 2p sp to 10p sp .) 
The cutoff prevents over-weighting of individual Monte Carlo elements so as to control the 
statistical error. 



At low temperatures, a very high density of small bubbles can coexist with a moderate density 
of large, low-pressure voids. Such distributions are most efficiently treated by making p max 
size-dependent, so that the maximum macroparticle densities are high in the region of bubbles, 
but low in the region of voids. Macroparticles can freely wander between the two regions. 
Accordingly, if macroparticle A moves to a region where pa > Pmax, it may be split into two 
identical parts; or if two MC elements at the same coordinate have Pa + Pb < Pmax, they may 
be united into one. 



In problems of reversible nucleation and growth, small MC clusters may shrink and disappear. 
It is computationally inefficient to follow unstable clusters by Monte Carlo methods. Accord- 
ingly, macroparticles of the smallest vacancy clusters (with both n vac < N™ n < N™ E and 
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nhei — 0) are deleted at the end of each time-step and their density returned to the correspond- 
ing element of Pme- (The numerical solution of the ME automatically accommodates any sub- 
sequent transients by adjusting its internal time-steps.) The minimum MC size should be large 
enough that macroparticles at the threshold only rarely shrink to monomer sizes during the 
interval r. It should also be far enough from that the cycle ME— >MC— >ME (involving 

creation of a new macroparticle, shrinkage of the constituent clusters, and transfer of that el- 
ement back to p ME ) is infrequent. In practice, N mm = N M£ /4 is used, and these two criteria 
are accomodated by taking the largest possible N ME . Helium clusters are never returned from 
MC to ME distribution; helium emission is not permitted, so the clusters will only grow along 
the helium axis. 

In the examples considered here, all ME clusters are sub-critical for ~ 60, so that newly- 
created MC particles frequently shrink and are annihilated. This is especially true at low 
temperatures, when the proliferation of small voids favors vacancy/interstitial recombination. 
Here, this "rare event problem" for nucleation of stable voids from small vacancy clusters is 
at least improved from conventional kinetic Monte Carlo methods, where even the monomers 
would be treated stochastically. Ultimately, direct application shows this mixed scheme is suit- 
able for radiation damage calculations to high doses. 



4.4. MC-MC reactions 



Coalescence problems are frequently treated by a Markov Monte Carlo method [15J . A straight- 
forward approach defines a finite volume, V, containing N (i.e., N MC ) discrete clusters of sizes 
{n} that stochastically evolve to a new N — 1 population {n'} through the binary coalescence 
of any pair of particles. The average rate of reaction between the zth and jth particles is simply 
A^(rij, iij)/V 2 per unit volume. The total rate of reaction of all N clusters is Rn, where 

i 

Ri = Yl R k,N (14) 

k=l 

and 

^ = £1%^ (15) 
k=i z 

in terms of the sum over reactions in the entire volume, V, assuming they are uncorrelated and 
occur in parallel. R^jy is proportional to the rate at which cluster i reacts with all other clusters. 

A stochastic sequence of discrete reactions may be constructed from these parameters. The 
random interval, £, to the next reaction is obtained from a uniform variate, x G (0, 1], as [49J: 

C = -\n(x)/R N (16) 

The first cluster of the random reaction pair, i, is selected with a probability proportional to 
Ri,N, from y G (0, 1] and 

Rj— i Rj , _ x 

ir < v - 7T (17) 

where Rq = 0. Finally, the reaction counterpart, j, is identified from z G (0, 1] and 
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with R it0 = 0. This selects j with a probability proportional to |i^(n i; n 3 -)/V. The procedure 
repeatedly selects new x, y, and 2 for the next event, increments the system time by £, performs 
the reaction i+j, and recalculates i? for the next iteration. This repeats until the elapsed time 
exceeds r. Since the last reaction falls outside the desired interval, it is discarded without being 
performed. The procedure may then be started anew for the next time-step. 

The choice of two random numbers to select the pair, i, j, differs from the usual scheme, where 
the pair is selected from a single value. In either case, the search for i and j takes o(log2{N)) 
operations using the method of bisection [50]. However, separate selection of i and j makes it 
possible to record all R m with o(N) storage space and a one-time computational effort of o(N 2 ). 
Once i is determined, i?i, m may be tabulated with o(N) effort for all m, so the full matrix need 
not be stored. Finally, after i and j react, the R m may be updated with o(N) effort by replacing 
only those terms involving the old clusters i and j with the results for a single new, coalesced 
cluster, and re-indexing to account for the lost cluster. Since Rn is an extensive quantity for 
a given total density, evolution of iV particles over a finite interval requires o(N 2 ) effort and 
o(N) storage. Specifying the binary reaction rate coefficients, K, as a half-triangular matrix 
increases the efficiency marginally. 

This MC scheme has difficulty modeling widely varying concentrations of reactants (e.g., the 
monomer density is typically orders of magnitude higher than the large clusters for radiation 
damage). Also, N decreases after every coalescence, which increases the statistical noise. There 
are methods that preserve N [T8f26] , but it is possible to encompass a wider range of densities at 
the same time. In the approach taken here, the discrete MC elements are macroparticles, widely 
used in, e.g., non-equilibrium simulations of plasma physics [IE]- (This is distinct from related, 
"weighted particle" schemes for coagulation []!?J22|) 51|.) Here, the j-th macroparticle in the sys- 
tem consists of an ensemble of clusters all of the same composition (n,, pj). Consistent with the 
Gillespie procedure, macroparticle reactions are evaluated discretely, so clusters in an ensem- 
ble react simultaneously but otherwise stochastically. However, here reactants will generally 
have different ensemble densities, pi < pn, which are independent of their sizes/compositions, 
tul, uh- The lower-density macroparticle, L, reacts completely, leaving behind an unchanged 
portion p H — p L of clusters from the higher-density ensemble, H. The total cluster density de- 
clines, but N stays constant, and iV-dependent errors remain steady over time. 



Macroparticle reaction rates (analogous to Eq. 15 ) are defined so as to reproduce the continuum 



limit as iV — > 00. Pairs i and j, with pi < pj, react according to: 

(Pi, Pi) + ( n i, Pj) ( n i + n j> Pi) + ( n j> Pj - Pi) (19) 

at an average rate of K(rii, tij)pj. Two macroparticles of the same density (pi = pj = p; i ^ j) 
react as: 

(n 4 , p) + (n„ p) -> (m + n„ p/2) + (m + n j; p/2) (20) 

at an average rate of K(rii, vij)p. The product is simply split into two equal pieces so that iV 
remains constant. Finally, the individual clusters within a single macroparticle ensemble may 
coalesce with one another, so there is also a unary reaction process: 

(m, Pi ) -> (2m, pi/2) (21) 
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which also proceeds at an average rate of K(rii, n^p^. This possibility modifies Eq. 15 to include 
a non-zero reaction rate for % — j. 



Macroparticle dynamics never corresponds to an atomistic simulation for finite N. Instead, this 
is a distinct, approximative discretization of the continuum equations themselves, in the same 
spirit as earlier approaches [5]. Again, P(t) is approximated here by a sparse set of elements 
without arbitrarily imposing some coarse-graining of finite difference equations for the distri- 
bution. Since the computational cost scales as o(N 2 ) for a dense reaction matrix, the method is 
also efficient. This is especially advantageous in higher dimensions, e.g., in describing helium- 
vacancy-impurity clusters. 



4.5. ME-MC reactions 



Additional schemes are required for treating reactions between ME and MC elements. In the 
continuum approximation, reaction with external entities, n ^ ME, introduces unary sink 
terms to the rate equation for m e ME, cf . Eqs. |9l [TO] 



S m (t)p m (t) 



Pm (t)U(N ME /2 - m) (22) 



J2 K(m,n(t))p n (t) + K(m,d)p d (t) 

LneAfC 

where the summation includes all elements {(n(i), p n {t)} MC at time t and where K(m, d) in- 
cludes reactions with network dislocations. The sink term, S, is identically zero for constrained 
ME defects. At present, K(m, d) is only nonzero for m = (1, 0), (—1, 0) and for vacancy emis- 
sion K(0, d). 

The counterpart to Eq.|22]is expressed for n e MC in the macroparticle scheme by: 

(n,p n ) A/C ^(m + n,p n ) MC (23) 
as a discrete reaction with an average rate of K(m, n)p m . A stochastic sequence of reactions at 



these average rates will approach the continuum behavior of Eq. [22jin the limit N MC — > oo. A 
single reaction can change a macroparticle size, cross-section, and reaction rate substantially, 
if m is comparable in size to n. Accordingly, ME-MC reactions for such "small" MC clusters 
are included by the Markov Monte Carlo scheme described above, and the reaction parameters 
are immediately updated to reflect the change, before evaluating the next reaction. 

Reaction events are randomly performed from the N MC x N ME matrix of reaction rates, at 
overall rate Q. If the next event occurs within the desired interval, the ith MC element is 
selected as a reactant with probability Qi/Q, where: 

Qi = Y^K{n i ,m j )p ni . (24) 
j 

for reactive elements nx, e ME and: 

N MC 

Q= EQi (25) 

i=l 

The jth ME element is selected as a reactant with probability K(rii, m.j)p mj /Qi- Finally, the 
time index is updated, the reaction is performed, and Q is revised. This is analogous to the 
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Markov procedure for MCMC reactions, except that the reaction matrix is full-rectangular 
rather than half-triangular and that the rates are always proportional to the density of the ME 
reactant. 



As for the corresponding evolution of F ME , the instantaneous source/sink terms, Eq. 22 , change 
after each discrete reaction event in P MC ( possibly multiple times during the interval r. It 
is not computationally practical to evolve Pme over each individual Markov sub-step, £, to 
account for this. Instead, Pme is evolved over the full time-step r by operator splitting, after 
all ME-MC and MC-MC reactions in Pmc are performed. To minimize any convergence error, 
the instantaneous sink strength can be replaced with a weighted time average over the interval: 



of 

On 



ast 
m 



t Q +T 



dt 



T Jt 
1 

r 



£ K(m, n(t))p m (t) 

neMC* 



^i [ :(m,n(t j _ 1 ))p m (t j _ 1 ) 



(26) 
(27) 



finally expressed as a sum over sub-intervals, between discrete reaction times, tj. 



Such attention to detail is unnecessary for large MC clusters (and for network dislocations), 
where rapid reactions with highly mobile defects (i.e., small m) do not substantially change 
the sink strength over short intervals. Thus, it is sufficient to update parameters for the large 
n clusters at the end of each time-step. In this case, MC clusters are evolved using a Poisson- 
distributed random variate, P(x), [20)152] for the number of reactions that occurs during r. 
These MC elements are only updated at to + T , with all reactions accumulated in each of the 
N ME channels: 



(n,Pr 



n 



Pr 



(2f 



mP TK(m.,n)p r 

Equation|28]is the discrete analogue of the Gaussian-distributed random walk used previously 
[5] . The corresponding ME sink term is: 



S s Jr(t )= ]T K(m,n(t ))p n (t ) + K( m ,d)p d (t ) 



(29) 



neMC 



including the dislocation contribution, assuming that Pd(t) is slowly changing. 



Finally, discrete reactions could also be evaluated by a rejection method, given a Majorant 
kernel M(m, n) > K(m, n) [51]. For example, the reaction rates, M, can be evaluated on a 
coarse grid of rij and all reactants < n < + 1 be treated alike. In another approach, M 
may be chosen to be a sum of products 



M(m,n)= M(m) ■ M(n) 



(30) 



It is then only required to evaluate N MC vectors, Ai, (of one or more dimensions) and to take 
dot products. Either approach is easier than directly computing N M (N MC + 1)/2 binary rate 
coefficients for Eqs. 14 , 15 The Majorant kernel is selected to be easy to evaluate and to predict 
a faster (or equal) event rate than the true system. To correct for any overestimate, the time 
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index is updated according to the usual Markov Monte Carlo procedure, but the reaction is 
only performed if a uniform variate, w G (0, 1] also satisfies w < K(m, n)/M(m, n). Thus, 
excess events predicted by M are rejected (with the required probability 1 — K/M). At present, 
the full reaction rate coefficients from Eq. [T] can be evaluated very efficiently, so this method 
is not employed here. However, it is expected to be advantageous when biased cavity-cavity, 
cavity-loop, and loop-loop reactions are included in the future. 



5. Results 



5.1. Monomer aggregation model 



A high density of trapping/recombination centers is believed to delay the onset of void swelling 
[%3|53|54j . Traps hinder void diffusion and coalescence and prolong the incubation stage. The 
simplest trapping model assumes that all dimers and larger clusters are immobile: D n = for all 
n £ {v , v h, h, i}, so that the last two summations in Eq. [3]are zero. If Eq. [2]is solved separately 
from the remainder of the Master Equation (Eq. [3) in a quasi-stationary approximation, then 
that problem may be solved by existing methods 



5|55] . However, here the problem is simply 



treated as a limit case of Smoluchowski's coagulation equation. 



Initial cluster populations are shown in Figs. [T]j3]for type-316 stainless steel irradiated to low 
doses at 10~ 6 dpa/s and 300, 500, and 700 C. It is well-known that helium-vacancy clusters 
may be separated into distinct species (of equilibrium bubbles and stable or unstable voids), 
according to their size-dependent free energies. Accordingly, the figures are marked with black 
lines where the net average vacancy addition rate for the defect clusters approaches zero. The 
leftmost black line in Figs. [T]j3] represents a hard wall for over-pressurized bubbles: by fiat, 
bubbles cannot reach densities greater than 2 helium per vacant site. Here, this is imposed by 
disallowing further reactions with helium- and self-interstitials. Other lines separate clusters 
that add or lose vacancies on average. Small, over-pressurized bubbles tend to add vacancies 
until reaching the next line in the Figures, where stable helium bubbles are in dynamic equilib- 
rium with the vacancy and interstitial population. (This approximates the line of bubbles with 
P ~ 7/2r, which would be in equilibrium in the absence of a vacancy and interstitial super- 
saturation.) The stability of these bubbles is reflected by their elevated concentration in that 
region, especially visible in Fig. [3] Stable bubbles tend to grow along the equilibrium line as 
they accumulate helium, while adjusting their vacancy content on average to remain in equi- 
librium. Finally, bubbles cannot exceed some critical helium content - larger clusters are stable 
voids that tend to add vacancies and grow to arbitrary size. This is seen in Fig. [2] there the 
clusters grow along the line of stable bubbles until reaching a critical helium content (11 he- 
liums), at which point they grow by adding vacancies in excess of helium, forming a plume of 
rapidly-growing voids in the size distribution. 

Voids are here simply taken to be cavities with higher vacancy /helium ratio than any bubble 
species of the same helium content. An approximately parabolic region under the black curves 
bounds a set of small, unstable voids that tend to lose vacancies and shrink back towards the 
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equilibrium bubble configuration. For example, this ranges from the origin to vacancy /helium 
compositions of (19,11) and (94,0) in Fig. [2} The rightmost solid line identifies the critical or 
unstable equilibrium void compositions; larger voids tend to add net vacancies with time. Note 
that a percentage of equilibrium bubbles in Fig.|2]are able to fluctuate in vacancy content across 
the barrier of unstable voids. That is, they become stable voids without having first reached 
the critical helium content. Similarly, helium dimers are readily able to cross the barrier of 
unstable voids in Fig. [TJ Very large voids ultimately become neutral (unbiased) sinks, adding 
helium/vacancies at an average rate of 1:200 (based on anticipated asymptotic swelling of 
1%/dpa and model helium generation around 50 appm/dpa). Thus, voids approach a line of 
constant composition. 

Except for a brief transient at the onset of irradiation, the vacancy monomer concentration 
decreases monotonically with time as the total sink strength of the microstructure rises with 
dose. After a few dpa, production of a-particles also peaks, and the helium monomer concen- 
tration also declines. During this extended period, equilibrium bubbles continue to grow by 
adding helium, they continue to reach the critical size, and they continue to become voids. 
However, the critical size for equilibrium bubbles increases with time (as a function of declining 
p v ), and the rate of formation of new helium dimer nuclei and bubble growth rates decrease (as 
a function of declining ph + Phv)- This causes the rate of void formation to decrease gradually 
with time, giving a broad void size distribution. Eventually, the larger stable bubbles become 
T EM- visible, and the overall size distribution becomes bimodal. 

The time-dependent volumetric swelling for this model is shown at a series of temperatures 
in Fig. |4j The low temperature system is initially dominated by large numbers of transient, 
unstable vacancy clusters (Fig.[T| that serve as recombination centers and suppress swelling. So 
many defect centers form that helium/ vacancy ratios are kept low, and helium plays a reduced 
role in the initial evolution. As a result, the visible cavity density (r > 0.5 nm) is sensitive to 
the surface energy parameter, 7: p V i S = 5 x 10 23 m -3 for 7(T) = 0.87o(T) and 1 x 10 24 m -3 for 
7(T) = O.570 (T). Eventually, some vacancy clusters acquire significant amounts of helium, and 
the system is filled with a high concentration of small equilibrium bubbles. These function as 
recombination centers; they keep the vacancy supersaturation low so that few, if any, bubbles 
grow into stable voids. They also keep the asymptotic swelling rate small. At and above 500 C, 
swelling is more a matter of helium bubble formation and growth towards critical sizes (Figs. [2] 
and|3|. The cavity density and swelling rates are therefore insensitive to 7. The steady swelling 
rate of 0.8 %/dpa at 500 C is consistent with void swelling in austenitic stainless steel [7f8] . 
At higher temperatures, the increased helium mobility results in fewer cavities (7-8 x 10 20 m~ 3 
at 700 C), and a smaller density of bubbles escape to become stable voids and contribute to 
steady swelling. 



5.2. Cluster coalescence model 



The other simplification of defect trapping is to neglect it entirely and assume that clusters 
diffuse freely according to their size. The predicted void size distribution changes significantly 
when coalescence is included. This is seen in Figs. [5] and |6j for the same temperatures as in 
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Figs. [T] and [2| Coalescence reactions continually, preferentially consume the smaller, more mo- 
bile clusters. The largest voids grow an order of magnitude larger through coalescence, making 
the distribution of stable void sizes substantially broader than before. Very large voids achieve 
such low diffusivities as to be effectively immobile; this results again in a terminal void popula- 
tion. At low temperatures, the removal of small unstable or equilibrium defect clusters reduces 
the number of recombination centers, suppresses damage annihilation, and speeds the forma- 
tion of large, stable voids. This enhances low temperature swelling. At high temperatures, this 
same coalescence of small clusters greatly reduces the total number of helium bubbles and voids, 
so that the total void sink strength is kept small and the asymptotic swelling rate is diminished 
compared to the monomer aggregation model (Fig. [7]). Small clusters are absorbed as rapidly 
as new ones form, which prevents the formation of a bimodal distribution of small equilibrium 
bubbles and large voids. These differences suggest that competition between trapping and coa- 
lescence of very small (mostly TEM-invisible) clusters significantly shapes the microstructure 
in real irradiated materials. 

When coalescence is included, the terminal void density and swelling rate remain sensitive to 
7 up to 500 C. The predicted void density at this temperature increases from 7 x 10 19 m -3 for 
7 = 0.757o to 7 x 10 20 m~ 3 for 7 = O.570. The swelling rate for the former case is only 0.3%/dpa 
at 50 dpa but reaches 0.8%/dpa for the latter. This suggests that either the cavity surface 
energy is substantially smaller than the value for the clean metal or that the vacancy clusters 
have much smaller mobilities than are modeled here. The swelling behavior finally becomes 
insensitive to the surface energy by 700C. In this limit, coalescence reduces the terminal void 
density to 4-5 x 10 18 m~ 3 . 



6. Conclusions 



This paper introduces a mixed Master Equation/Monte Carlo treatment of rate theory calcu- 
lations in a mean field continuum approximation. This enables flexible treatment of the defect 
density variables, using different algorithms to treat the various reactions as efficiently as possi- 
ble. The approximately quasi- stationary distribution of small, unstable or transient clusters is 
treated (as much as possible) by solving continuum rate equations. This eliminates the need to 
evaluate rapid individual reactions that mostly cancel one another. Larger clusters are treated 
by Monte Carlo methods, which treats clusters of arbitrary size and composition without the 
need for a fixed grid or artificial discretization of the defect distribution. A Markov method 
for smaller clusters accurately simulates rapid fluctuations in size and in the reaction parame- 
ters, and a Poisson-distributed random walk efficiently treats the more gradual evolution of the 
largest clusters. Finally, a macroparticle approach is introduced to encompass large differences 
in species densities in the Monte Carlo distribution. 

This hybrid scheme readily treats void/bubble evolution to high cumulative fluxes for temper- 
atures and dose rates that are characteristic of real reactor systems. Calculations demonstrate 
that void coalescence provides an important channel for consolidating vacancy defects into 
large, stable voids, controlling the duration of incubation and the terminal void density. It is ex- 
pected that thermal and radiation-induced micro-chemical evolution of solute and precipitate 
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distributions will influence the cluster mobility and thereby the macroscopic incubation and 
steady-swelling behavior. Some degree of void/bubble trapping seems to be required in order to 
obtain a bimodal bubble/void size distribution, while some coalescence may be needed to give a 
realistically low terminal void density at higher temperature. The cavity surface energy deter- 
mines the barrier for nucleation of stable voids and hence also affects the incubation behavior; 
this contribution becomes temperature- and time-dependent if oxygen is explicitly modeled. 
All of these effects can be addressed, in principle, by extensions of the method described here. 

These calculations also suggest the importance of additional, competing processes that are not 
evaluated at present, such as interstitial-interstitial aggregation or cluster annihilation from 
void-dislocation reaction. The methods described here can be extended to treat coalescence 
of loops as easily as voids, given a suitable binary reaction kernel. Such reactions should be 
included for reasons of consistency, besides their likely contribution to transient and steady 
swelling behavior. They would be especially important if radiation damage were introduced as 
a variety of pre-formed defect clusters. Based on the preliminary findings for cavity coalescence, 
more general defect cluster reactions are expected to have a significant influence on radiation 
swelling behavior. 
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Fig. 1. A portion of the void/bubble distribution for a model with mobile monomer defects and sessile clusters, 7 = 0.87o(T), 
at T=300 C, 10 — 6 dpa/s, and 32 X 10 — 3 dpa. The largest void in the distribution contains 110 vacancies. The solid lines 
display the loci of stable and unstable equilibrium cluster compositions, based on average vacancy accumulation rates. This 
distribution has not been smoothed - the pixellated appearance reflects discrete cluster compositions. (Fig2) 




# Vacancies 

Fig. 2. A portion of the void/bubble distribution as in Fig.[l] but at T=500 C, 1CT 6 dpa/s, and 16 X 10~ 3 dpa. The distribution 
has been smoothed for the large clusters, where Monte Carlo data is increasingly sparse. The solid lines display the stable and 
unstable equilibrium cluster compositions. 
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Fig. 3. The full void/bubble distribution as in Fig. [5] but at T=700 C. The curved solid line locates the stable equilibrium 
bubbles; the critical void size is not visible on this scale. 
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Fig. 4. Volumetric swelling curves versus dose in the model that excludes void-void coalescence. The cavity surface energy is 
fixed at 7(T) = 0.47 (T) 
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Fig. 5. A portion of the void/bubble distribution as in Fig. ^(300 C), but including void coalescence and with 7(T) = 0.57o(T). 
The distribution has been smoothed for the large clusters, where Monte Carlo data is sparse. 
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Fig. 6. The full void/bubble distribution as in Fig. [3] (500 C), but including void-void coalescence and with j(T) = 0.57o(T). 
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Fig. 7. Volumetric swelling curves versus dose in the model that includes void-void coalescence. The cavity surface energy is 
set to 7(T) = 0.47 (T). 
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Table 1 

Model material parameters for typc-316 stainless steel. 



Bulk parameters: 


Lattice constant ao 
Burgers vector b 
Atomic volume f! 
Shear modulus fi 
Poisson's ratio v 
Cascade efficiency ^Frenkel 


3.639 x 10- 10 m 

ao / V £ 

ao 3 /4 
8.295 x 10 10 Pa 

0.264 
n 1 

U. L 


Vacancy parameters: 


Relaxation volume 
Migration energy E m 
Formation energy Ej 
j? kjl maiioii ciiLiopy Of 
Pre-exponential factor 

Shear polarizability 


-0.2 n 

2.08 x 10~ 19 J 
2.88 x 10~ 19 J 

1.29 X 10~ 6 m 2 /s 
-2.4 X 10~ 18 


Self-interstitial 
parameters: 


Relaxation volume 

l\/1 1 <TT" Ct t" 1 rwi nY\t~*v<T\r fi 

iviigi aiiion ciieigy f-^m 
Pre-exponential factor 
Shear polarizability 


1.50 n 

^90 v 1 — 19 T 
1.29 X 10~ 6 m 2 /s 
-2.535 x 10~ 17 


Interstitial helium 
parameters: 


Relaxation volume 
Migration energy E m 
Pre-exponential factor 

oiicdii jjuidi izjciuiii l y 


0.60 n 
0.320 x 10~ 19 J 
1.29 x 10~ 6 m 2 /s 
—9 v 1 n — 17 


Substitutional helium 
parameters: 


Relaxation volume 
Migration energy E m 
Pre-exponential factor 

Shear polarizability 


-0.2 n 

2.08 x 10~ 19 J 
1.29 x 10~ 6 m 2 /s 
-2.4 x 10~ 18 


Void parameters: 


Relaxation volume 
Surface energy 70 (T = 0) 
Temperature derivative dyo/dT 
Adsorption factor A 
Migration energy E m 
Pre-exponential factor 




2.408 J/m 2 
0.440 X 10~ 3 J/m 2 /K 
0.45-0.80 
2.08 x 10~ 19 J 
1.29 X 10- 6 m 2 /s /ni /S 


Environmental 
parameters: 


Temperature T 
Flux <j> 
Damage efficiency £ 


300-700 C 
10~ 6 dpa/s 
0.1 
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